﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

using System.Windows.Forms;

using OSGeo.GDAL;


namespace HATerrain
{

    /// <summary>
    /// 用于表示地形数据的读取部分
    /// </summary>
    public  class TerrainFileHelper
    {

        public static TerrainGrid ReadHgtFile(string filename)
        {


            TerrainGrid terrain = new TerrainGrid();



            string _fileName;
             OSGeo.GDAL.Driver Driver;
         
             OSGeo.GDAL.Dataset Dataset;
         
             OSGeo.GDAL.DataType DataType;
        
             OSGeo.GDAL.Band Band;


            OSGeo.GDAL.Gdal.AllRegister();

            _fileName = filename;

            if (System.IO.File.Exists(filename) == false)
            {
                return null;
            }

            Dataset = OSGeo.GDAL.Gdal.Open(_fileName, Access.GA_ReadOnly);

            Driver = Dataset.GetDriver();

            Band = Dataset.GetRasterBand(1);

            DataType = Band.DataType;


            ///开始读取关键信息


            terrain.Info .SizeX = Dataset.RasterXSize;
            terrain.Info.SizeY = Dataset.RasterYSize;
            terrain.Info.Projection = Dataset.GetProjection();

            Dataset.GetGeoTransform(terrain.Info.Affine);    //映射

            //ASC文件中 CellSizeX与CellSizeY的绝对值是相等的
            terrain.Info.CellSizeX = Math.Abs(terrain.Info.Affine[1]);
            terrain.Info.CellSizeY = Math.Abs(terrain.Info.Affine[5]);

            //西北角坐标
            terrain.Info.CoordinateNW.X = terrain.Info.Affine[0];
            terrain.Info.CoordinateNW.Y = terrain.Info.Affine[3]; // this.CoordinateNE.Y + this.SizeY * Math.Abs(this.CellSizeY);

            //西南角坐标
            terrain.Info.CoordinateSW.X = terrain.Info.Affine[0];
            terrain.Info.CoordinateSW.Y = terrain.Info.CoordinateNW.Y - (terrain.Info.SizeY - 1) * Math.Abs(terrain.Info.CellSizeY);

            //东北角坐标
            terrain.Info.CoordinateNE.X = terrain.Info.CoordinateNW.X + (terrain.Info.SizeX - 1) * Math.Abs(terrain.Info.CellSizeX);
            terrain.Info.CoordinateNE.Y = terrain.Info.CoordinateNW.Y;
            //东南角坐标
            terrain.Info.CoordinateSE.X = terrain.Info.CoordinateNE.X;
            terrain.Info.CoordinateSE.Y = terrain.Info.CoordinateSW.Y;

            int hasval = 0;

            //最大、最小、无效值
            double max, min, nodata;

            double mean, std;

            Band.GetStatistics(0, 1, out min, out max, out mean, out std);

            //下面这两个怎么不正确，还是是别的意思
            //base.Band.GetMaximum(out max, out hasval);
            terrain.Info.DataMax = max;

            //base.Band.GetMinimum(out min, out hasval);
            terrain.Info.DataMin = min;

            //无效值

            

            Band.GetNoDataValue(out nodata, out hasval);
            terrain.Info.NoDataValue = nodata;



            //读取数据,首先判断数据类型
            dynamic _data = null;
            switch (Band .DataType)
            {
                case DataType.GDT_Byte:
                    _data = new byte[terrain.Info.SizeX * terrain.Info.SizeY];
                    break;
                case DataType.GDT_CFloat32:
                    _data = new float[terrain.Info.SizeX * terrain.Info.SizeY];
                    break;
                case DataType.GDT_CFloat64:
                    _data = new float[terrain.Info.SizeX * terrain.Info.SizeY];
                    break;
                case DataType.GDT_CInt16:
                    _data = new Int16[terrain.Info.SizeX * terrain.Info.SizeY];
                    break;
                case DataType.GDT_CInt32:
                    _data = new Int32[terrain.Info.SizeX * terrain.Info.SizeY];
                    break;
                case DataType.GDT_Float32:
                    _data = new float[terrain.Info.SizeX * terrain.Info.SizeY];
                    break;
                case DataType.GDT_Float64:
                    _data = new float[terrain.Info.SizeX * terrain.Info.SizeY];
                    break;
                case DataType.GDT_Int16:
                    _data = new Int16[terrain.Info.SizeX * terrain.Info.SizeY];
                    break;
                case DataType.GDT_Int32:
                    _data = new Int32[terrain.Info.SizeX * terrain.Info.SizeY];
                    break;
                case DataType.GDT_TypeCount:
                    break;
                case DataType.GDT_UInt16:
                    _data = new UInt16[terrain.Info.SizeX * terrain.Info.SizeY];
                    break;
                case DataType.GDT_UInt32:
                    _data = new UInt32[terrain.Info.SizeX * terrain.Info.SizeY];
                    break;
                case DataType.GDT_Unknown:
                    break;
                default:
                    break;
            }

             


            


             //因为对Band还没有深入的了解，先这样写，下来再研究改正
             if (Band != null&&_data!=null)
             {
                 Band.ReadRaster(0,0,terrain.Info.SizeX, terrain.Info.SizeY, _data, terrain.Info.SizeX, terrain.Info.SizeY, terrain.Info.PixelSpace , terrain.Info. LineSpace);
             }
             else
             {
                 //这个是用Dataset读栅格数据，下来再加并测试
                 //base .Dataset .ReadRaster ( xOff,  yOff,  sizeX,  sizeY,this .data ,buf_sizeX ,buf_sizeY ,1,
             }

             Dataset.FlushCache();

        
            //////////////////////////////////一下代码需考虑_data未被赋值的情况
           
            //TD.Data = new float[this.Dataset.RasterYSize * this.Dataset.RasterXSize];

            //for (int i = 0; i < this._data.Length; i++)
            //{
            //    TD.Data[i] = (float)this._data[i];
            //}

            //TD.Data2D = new float[TD.Header.SizeY][];

            //for (int i = 0; i < TD.Header.SizeY; i++)
            //{
            //    TD.Data2D[i] = new float[TD.Header.SizeX];
            //    for (int j = 0; j < TD.Header.SizeX; j++)
            //    {
            //        TD.Data2D[i][j] = TD.Data[i * TD.Header.SizeX + j];
            //    }
            //}

            //this._data = null;







             return terrain;


        }




        # region 


        
        public static TerrainData Open(string filename)
        {
            if (System.IO.File.Exists(filename) == false)
            {
                return null;
            }

            TerrainData TD = new TerrainData();

            string _fileName;
            OSGeo.GDAL.Driver Driver;

            OSGeo.GDAL.Dataset Dataset;

            OSGeo.GDAL.DataType DataType;

            OSGeo.GDAL.Band Band;


            OSGeo.GDAL.Gdal.AllRegister();

            _fileName = filename;

            if (System.IO.File.Exists(filename) == false)
            {
                return null;
            }

            Dataset = OSGeo.GDAL.Gdal.Open(_fileName, Access.GA_ReadOnly);

            Driver = Dataset.GetDriver();

            Band = Dataset.GetRasterBand(1);

            DataType = Band.DataType;

            TD.Header = ReadHeader(Dataset);


            //读取数据,首先判断数据类型
            dynamic _data = null;
            switch (Band.DataType)
            {
                case DataType.GDT_Byte:
                    _data = new byte[TD.Header.SizeX * TD.Header.SizeY];
                    TD.Data2D = new byte[TD.Header.SizeY, TD.Header.SizeX];
                    break;
                case DataType.GDT_CFloat32:
                    _data = new float[TD.Header.SizeX * TD.Header.SizeY];
                    TD.Data2D = new float[TD.Header.SizeY, TD.Header.SizeX];
                    break;
                case DataType.GDT_CFloat64:
                    _data = new float[TD.Header.SizeX * TD.Header.SizeY];
                    TD.Data2D = new float[TD.Header.SizeY, TD.Header.SizeX];
                    break;
                case DataType.GDT_CInt16:
                    _data = new Int16[TD.Header.SizeX * TD.Header.SizeY];
                    TD.Data2D = new Int16[TD.Header.SizeY, TD.Header.SizeX];
                    break;
                case DataType.GDT_CInt32:
                    _data = new Int32[TD.Header.SizeX * TD.Header.SizeY];
                    TD.Data2D = new Int32[TD.Header.SizeY, TD.Header.SizeX];
                    break;
                case DataType.GDT_Float32:
                    _data = new float[TD.Header.SizeX * TD.Header.SizeY];
                    TD.Data2D = new float[TD.Header.SizeY, TD.Header.SizeX];
                    break;
                case DataType.GDT_Float64:
                    _data = new float[TD.Header.SizeX * TD.Header.SizeY];
                    TD.Data2D = new float[TD.Header.SizeY, TD.Header.SizeX];
                    break;
                case DataType.GDT_Int16:
                    _data = new Int16[TD.Header.SizeX * TD.Header.SizeY];
                    TD.Data2D = new Int16[TD.Header.SizeY, TD.Header.SizeX];
                    break;
                case DataType.GDT_Int32:
                    _data = new Int32[TD.Header.SizeX * TD.Header.SizeY];
                    TD.Data2D = new Int32[TD.Header.SizeY, TD.Header.SizeX];
                    break;
                case DataType.GDT_TypeCount:
                    break;
                case DataType.GDT_UInt16:
                    _data = new UInt16[TD.Header.SizeX * TD.Header.SizeY];
                    TD.Data2D = new UInt16[TD.Header.SizeY, TD.Header.SizeX];
                    break;
                case DataType.GDT_UInt32:
                    _data = new UInt32[TD.Header.SizeX * TD.Header.SizeY];
                    TD.Data2D = new UInt32[TD.Header.SizeY, TD.Header.SizeX];
                    break;
                case DataType.GDT_Unknown:
                    break;
                default:
                    break;
            }
            


            //因为对Band还没有深入的了解，先这样写，下来再研究改正
            if (Band != null && _data != null)
            {
                Band.ReadRaster(0, 0, TD.Header.SizeX, TD.Header.SizeY, _data, TD.Header.SizeX, TD.Header.SizeY, TD.Header.PixelSpace, TD.Header.LineSpace);
            }
            else
            {
                //这个是用Dataset读栅格数据，下来再加并测试
                //base .Dataset .ReadRaster ( xOff,  yOff,  sizeX,  sizeY,this .data ,buf_sizeX ,buf_sizeY ,1,
            }

            Dataset.FlushCache();

            if (_data == null)  //这里怎么判断呢
            {
                TD.Data2D=null;
            }
            else
            {
                for (int i = 0; i < TD.Header.SizeY; i++)
                {
                    for (int j = 0; j < TD.Header.SizeX; j++)
                    {
                        TD.Data2D[i,j] = _data[i * TD.Header.SizeX + j];                        
                    }
                }
            }


            return TD;

        }


        /// <summary>
        /// 读取地形文件的头文件
        /// </summary>
        /// <returns></returns>
        public static TerrainHeader ReadHeader(OSGeo.GDAL.Dataset Dataset)
        {

            OSGeo.GDAL.Driver Driver = Dataset.GetDriver();

            OSGeo.GDAL.Band Band = Dataset.GetRasterBand(1);

            TerrainHeader TH = new TerrainHeader();


            TH.SizeX = Dataset.RasterXSize;
            TH.SizeY = Dataset.RasterYSize;
            TH.Projection = Dataset.GetProjection();

            Dataset.GetGeoTransform(TH.Affine);    //映射

            //ASC文件中 CellSizeX与CellSizeY的绝对值是相等的
            TH.CellSizeX = Math.Abs(TH.Affine[1]);
            TH.CellSizeY = Math.Abs(TH.Affine[5]);

            //西北角坐标
            TH.CoordinateNW.X = TH.Affine[0];
            TH.CoordinateNW.Y = TH.Affine[3]; // this.CoordinateNE.Y + this.SizeY * Math.Abs(this.CellSizeY);

            //西南角坐标
            TH.CoordinateSW.X = TH.Affine[0];
            TH.CoordinateSW.Y = TH.CoordinateNW.Y - (TH.SizeY - 1) * Math.Abs(TH.CellSizeY);

            //东北角坐标
            TH.CoordinateNE.X = TH.CoordinateNW.X + (TH.SizeX - 1) * Math.Abs(TH.CellSizeX);
            TH.CoordinateNE.Y = TH.CoordinateNW.Y;
            //东南角坐标
            TH.CoordinateSE.X = TH.CoordinateNE.X;
            TH.CoordinateSE.Y = TH.CoordinateSW.Y;

            int hasval = 0;

            //最大、最小、无效值
            double max, min, nodata;

            double mean, std;

            Band.GetStatistics(0, 1, out min, out max, out mean, out std);

            //下面这两个怎么不正确，还是是别的意思
            //base.Band.GetMaximum(out max, out hasval);
            TH.DataMax = max;

            //base.Band.GetMinimum(out min, out hasval);
            TH.DataMin = min;

            //无效值
            Band.GetNoDataValue(out nodata, out hasval);
            TH.NoDataValue = nodata;

            return TH;

        }


        /// <summary>
        /// 合并两个地形文件
        /// </summary>
        /// <param name="data1"></param>
        /// <param name="data2"></param>
        /// <returns></returns>
        public static TerrainData Union(TerrainData data1, TerrainData data2)
        {

            if (data1 == null || data1.Data2D == null || data1.Data2D.Length <= 0)
            {
                if (data2 == null || data2.Data2D == null || data2.Data2D.Length <= 0)
                {
                    return null;
                }
                else
                {
                    return data2;
                }
            }
            else
            {
                if (data2 == null || data2.Data2D == null || data2.Data2D.Length <= 0)
                {
                    return data1;
                }
                else
                {

                }
            }



            if (data1.Header.Projection != data2.Header.Projection)
            {
                //如果投影不相同，则不能合并，
                //return null;
            }

            //定义一个新的TifManager——合成后的TifManager
            TerrainData newData = new TerrainData();


            //先找出两个文件的最大最小坐标
            //用这个方法转换坐标，必须在debug目录下存在proj.dll这个控件！我用的是Mapwindow带的这个控件
            OSGeo.OSR.CoordinateTransformation CT = new OSGeo.OSR.CoordinateTransformation
                (new OSGeo.OSR.SpatialReference(data2.Header.Projection), new OSGeo.OSR.SpatialReference(data1.Header.Projection));

            double[] ptOut = new double[3];

            CT.TransformPoint(ptOut, data2.Header.CoordinateSW.X, data2.Header.CoordinateSW.Y, 0);

            //=====目前对于我们的应用来说，上面这几句代码(坐转的转换)基本没什么用====================================//

            double maxX, minX, maxY, minY;

            minX = Math.Min(data1.Header.CoordinateSW.X, ptOut[0]);
            minY = Math.Min(data1.Header.CoordinateSW.Y, ptOut[1]);

            CT.TransformPoint(ptOut, data2.Header.CoordinateNE.X, data2.Header.CoordinateNE.Y, 0);

            maxX = Math.Max(data1.Header.CoordinateNE.X, ptOut[0]);
            maxY = Math.Max(data1.Header.CoordinateNE.Y, ptOut[1]);


            newData.Header.CoordinateNW = new HAPoint(minX, maxY, 0);
            newData.Header.CoordinateNE = new HAPoint(maxX, maxY, 0);
            newData.Header.CoordinateSE = new HAPoint(maxX, minY, 0);
            newData.Header.CoordinateSW = new HAPoint(minX, minY, 0);


            //网格间距
            double CellSizeX = 0, CellSizeY = 0;
            CellSizeX = Math.Min(Math.Abs(data1.Header.CellSizeX), Math.Abs(data2.Header.CellSizeX));
            CellSizeY = Math.Min(Math.Abs(data1.Header.CellSizeY), Math.Abs(data2.Header.CellSizeY));

            newData.Header.CellSizeX = CellSizeX;
            newData.Header.CellSizeY = CellSizeY;

            //文件行列数
            int rows = 0, cols = 0;
            cols = (int)((maxX - minX) / CellSizeX) + 1;
            rows = (int)((maxY - minY) / CellSizeY) + 1;

            //if (cols < data1.Header.SizeX + data2.Header.SizeX - 1)
            //{
            //    cols = data1.Header.SizeX + data2.Header.SizeX;
            //}
            //if (rows < data1.Header.SizeY + data2.Header.SizeY - 1)
            //{
            //    rows = data1.Header.SizeY + data2.Header.SizeY;
            //}

            newData.Header.SizeX = cols;
            newData.Header.SizeY = rows;

            //投影,暂时还不会写，下来再研究
            newData.Header.Projection = data1.Header.Projection;

            //string str="";
            //int l = OSGeo.OSR.Osr.GetWellKnownGeogCSAsWKT("WGS84", out  str);

            //映射
            double[] Geo = new double[6];
            Geo[0] = minX;
            Geo[1] = CellSizeX;
            Geo[2] = 0;
            Geo[3] = maxY;
            Geo[4] = 0;
            Geo[5] = 0 - CellSizeY;

            newData.Header.Affine = Geo;

            //无效数据
            newData.Header.NoDataValue = data1.Header.NoDataValue;

            //数据
            //newData.Data = new double[rows * cols];
            int[] Data = new int[rows * cols];

            for (int i = 0; i < Data.Length; i++)
            {
                Data[i] = (int)newData.Header.NoDataValue;
            }

            //先给newData导入第一个地形的数据

            int startX1, startY1;
            startX1 = (int)((data1.Header.CoordinateNW.X - newData.Header.CoordinateNW.X) / CellSizeX);
            startY1 = (int)((newData.Header.CoordinateNW.Y - data1.Header.CoordinateNW.Y) / CellSizeY);

            for (int i = 0; i < data1.Header.SizeY; i++)
            {
                for (int j = 0; j < data1.Header.SizeX; j++)
                {
                    if ((startY1 + i) * newData.Header.SizeX + startX1 + j < Data.Length)
                    {
                        Data[(startY1 + i) * newData.Header.SizeX + startX1 + j] = data1.Data2D[i,j];
                    }

                    //测试
                    //CT.TransformPoint(ptOut, otherTif.CoordinateSW.X, otherTif.CoordinateSW.Y, 0);
                }
            }

            //给newTif导入第二个地形的数据

            int startX2, startY2;

            CT.TransformPoint(ptOut, data2.Header.CoordinateNW.X, data2.Header.CoordinateNW.Y, 0);


            startX2 = (int)((ptOut[0] - newData.Header.CoordinateSW.X) / CellSizeX);
            startY2 = (int)((newData.Header.CoordinateNW.Y - ptOut[1]) / CellSizeY);

            for (int i = 0; i < data2.Header.SizeY; i++)
            {
                for (int j = 0; j < data2.Header.SizeX; j++)
                {
                    if ((startY2 + i) * newData.Header.SizeX + startX2 + j < Data.Length)
                    {
                        Data[(startY2 + i) * newData.Header.SizeX + startX2 + j] = data2.Data2D[i,j];
                    }

                }
            }

            newData.Data2D = new dynamic[newData.Header.SizeY, newData.Header.SizeX];

            for (int i = 0; i < newData.Header.SizeY; i++)
            {
                for (int j = 0; j < newData.Header.SizeX; j++)
                {
                    newData.Data2D[i,j] = Data[i * newData.Header.SizeX + j];
                }
            }

            return newData;

        }

        /// <summary>
        /// 根据utm范围从经纬度文件裁剪
        /// </summary>
        /// <param name="TD"></param>
        /// <param name="minx"></param>
        /// <param name="miny"></param>
        /// <param name="maxx"></param>
        /// <param name="maxy"></param>
        /// <param name="utmZone"></param>
        /// <returns></returns>
        public static TerrainData Clip(TerrainData TD, double minx, double miny, double maxx, double maxy, int utmZone)
        {
            TerrainData newTD = new TerrainData();

            //新文件的四个点坐标
            newTD.Header.CoordinateNE = new HAPoint(maxx, maxy, 0);
            newTD.Header.CoordinateNW = new HAPoint(minx, maxy, 0);
            newTD.Header.CoordinateSE = new HAPoint(maxx, miny, 0);
            newTD.Header.CoordinateSW = new HAPoint(minx, miny, 0);

            //===================================================================================//

            OSGeo.OSR.SpatialReference utmPrj = new OSGeo.OSR.SpatialReference("");
            utmPrj.SetProjCS("UTM " + utmZone.ToString() + "(WGS84)");
            utmPrj.SetWellKnownGeogCS("WGS84");
            if (TD.Header.CoordinateSW.Y > 0)
            {
                utmPrj.SetUTM(utmZone, 1); //北半球
            }
            else
            {
                utmPrj.SetUTM(utmZone, 0); //南半球
            }
            //经纬度空间投影
            OSGeo.OSR.SpatialReference geoPrj = new OSGeo.OSR.SpatialReference(TD.Header.Projection);

            //坐标转换工具
            OSGeo.OSR.CoordinateTransformation utmToLonglat = new OSGeo.OSR.CoordinateTransformation(utmPrj, geoPrj);
            OSGeo.OSR.CoordinateTransformation longlatToUtm = new OSGeo.OSR.CoordinateTransformation(geoPrj, utmPrj);


            //新文件的投影
            string prj = "";
            utmPrj.ExportToWkt(out prj);
            newTD.Header.Projection = prj;


            //将原文件的经纬度坐标转换为utm坐标

            double utmx, utmy;
            double[] coorOut = new double[3];

            //经纬度转换为utm坐标
            //西南角坐标
            utmx = TD.Header.CoordinateSW.X;
            utmy = TD.Header.CoordinateSW.Y;

            longlatToUtm.TransformPoint(coorOut, utmx, utmy, 0);
            HAPoint ptSW = new HAPoint(coorOut[0], coorOut[1], 0);


            //西北角坐标
            utmx = TD.Header.CoordinateNW.X;
            utmy = TD.Header.CoordinateNW.Y;

            longlatToUtm.TransformPoint(coorOut, utmx, utmy, 0);
            HAPoint ptNW = new HAPoint(coorOut[0], coorOut[1], 0);


            //东北角坐标
            utmx = TD.Header.CoordinateNE.X;
            utmy = TD.Header.CoordinateNE.Y;

            longlatToUtm.TransformPoint(coorOut, utmx, utmy, 0);
            HAPoint ptNE = new HAPoint(coorOut[0], coorOut[1], 0);

            //东南角坐标
            utmx = TD.Header.CoordinateSE.X;
            utmy = TD.Header.CoordinateSE.Y;

            longlatToUtm.TransformPoint(coorOut, utmx, utmy, 0);
            HAPoint ptSE = new HAPoint(coorOut[0], coorOut[1], 0);

            double utmMinX = 0, utmMaxX = 0, utmMinY = 0, utmMaxY = 0;

            utmMinX = Math.Min(ptNW.X, ptSW.X);
            utmMaxX = Math.Max(ptNE.X, ptSE.X);
            utmMinY = Math.Min(ptSW.Y, ptSE.Y);
            utmMaxY = Math.Max(ptNW.Y, ptNE.Y);

            //===============原文件的坐标转换完毕==========================================

            //判断参数的裁切范围是否在文件的范围内
            if (minx < utmMinX || miny < utmMinY || maxx > utmMaxX || maxy > utmMaxY)
            {
                return null;
            }

            //新文件的网格间距
            double cellsizeX = (utmMaxX - utmMinX) / TD.Header.SizeX;
            double cellsizeY = (utmMaxY - utmMinY) / TD.Header.SizeY;

            newTD.Header.CellSizeX = cellsizeX;
            newTD.Header.CellSizeY = cellsizeY;

            //======================以上只是为了求出新文件的网格间距============================

            int sizeX = (int)Math.Round((maxx - minx) / cellsizeX);
            int sizeY = (int)Math.Round((maxy - miny) / cellsizeY);

            newTD.Header.SizeX = sizeX;
            newTD.Header.SizeY = sizeY;

            //====================以上是新文件的网格数=============================================

            //新文件的映射
            double[] geo = new double[6];
            geo[0] = minx;
            geo[1] = cellsizeX;
            geo[2] = 0;
            geo[3] = maxy;
            geo[4] = 0;
            geo[5] = 0 - Math.Abs(cellsizeY);

            newTD.Header.Affine = geo;

            //求新文件的高程数据
            int[] Data = new int[sizeX * sizeY];

            for (int j = 0; j < sizeY; j++)
            {
                for (int i = 0; i < sizeX; i++)
                {
                    double X = minx + i * cellsizeX;
                    double Y = miny + j * cellsizeY;
                    //转换成经纬度
                    utmToLonglat.TransformPoint(coorOut, X, Y, 0);

                    //插值算法，求出该点的高程值
                    Data[j * sizeX + i] = (int)GetPtElevation(TD, coorOut[0], coorOut[1]);
                }
            }

            newTD.Data2D = new dynamic[newTD.Header.SizeY, newTD.Header.SizeX];

            for (int i = 0; i < newTD.Header.SizeY; i++)
            {
                for (int j = 0; j < newTD.Header.SizeX; j++)
                {
                    newTD.Data2D[i,j] = Data[i * newTD.Header.SizeX + j];
                }
            }


            return newTD;
        }



        public static TerrainData Clip(TerrainData TD, double minx, double miny, double maxx, double maxy, int Lat1, int Lat2, int Lat, int Long,double falseEasting,double falseNorthing)
        {
            TerrainData newTD = new TerrainData();

            //新文件的四个点坐标
            newTD.Header.CoordinateNE = new HAPoint(maxx, maxy, 0);
            newTD.Header.CoordinateNW = new HAPoint(minx, maxy, 0);
            newTD.Header.CoordinateSE = new HAPoint(maxx, miny, 0);
            newTD.Header.CoordinateSW = new HAPoint(minx, miny, 0);

            //===================================================================================//

            OSGeo.OSR.SpatialReference lccPrj = new OSGeo.OSR.SpatialReference("");
            lccPrj.SetProjCS("LCC" + "(WGS84)");
            lccPrj.SetWellKnownGeogCS("WGS84");

            lccPrj.SetLCC(Lat1, Lat2, Lat, Long, falseEasting, falseNorthing);
          
            //经纬度空间投影
            OSGeo.OSR.SpatialReference geoPrj = new OSGeo.OSR.SpatialReference(TD.Header.Projection);

            //坐标转换工具
            OSGeo.OSR.CoordinateTransformation lccToLonglat = new OSGeo.OSR.CoordinateTransformation(lccPrj, geoPrj);
            OSGeo.OSR.CoordinateTransformation longlatToLcc= new OSGeo.OSR.CoordinateTransformation(geoPrj, lccPrj);


            //新文件的投影
            string prj = "";
            lccPrj.ExportToWkt(out prj);
            newTD.Header.Projection = prj;


            //将原文件的经纬度坐标转换为utm坐标

            double utmx, utmy;
            double[] coorOut = new double[3];

            //经纬度转换为utm坐标
            //西南角坐标
            utmx = TD.Header.CoordinateSW.X;
            utmy = TD.Header.CoordinateSW.Y;

            longlatToLcc.TransformPoint(coorOut, utmx, utmy, 0);
            HAPoint ptSW = new HAPoint(coorOut[0], coorOut[1], 0);


            //西北角坐标
            utmx = TD.Header.CoordinateNW.X;
            utmy = TD.Header.CoordinateNW.Y;

            longlatToLcc.TransformPoint(coorOut, utmx, utmy, 0);
            HAPoint ptNW = new HAPoint(coorOut[0], coorOut[1], 0);


            //东北角坐标
            utmx = TD.Header.CoordinateNE.X;
            utmy = TD.Header.CoordinateNE.Y;

            longlatToLcc.TransformPoint(coorOut, utmx, utmy, 0);
            HAPoint ptNE = new HAPoint(coorOut[0], coorOut[1], 0);

            //东南角坐标
            utmx = TD.Header.CoordinateSE.X;
            utmy = TD.Header.CoordinateSE.Y;

            longlatToLcc.TransformPoint(coorOut, utmx, utmy, 0);
            HAPoint ptSE = new HAPoint(coorOut[0], coorOut[1], 0);

            double utmMinX = 0, utmMaxX = 0, utmMinY = 0, utmMaxY = 0;

            utmMinX = Math.Min(ptNW.X, ptSW.X);
            utmMaxX = Math.Max(ptNE.X, ptSE.X);
            utmMinY = Math.Min(ptSW.Y, ptSE.Y);
            utmMaxY = Math.Max(ptNW.Y, ptNE.Y);

            //===============原文件的坐标转换完毕==========================================

            //判断参数的裁切范围是否在文件的范围内
            if (minx < utmMinX || miny < utmMinY || maxx > utmMaxX || maxy > utmMaxY)
            {
                return null;
            }

            //新文件的网格间距
            double cellsizeX = (utmMaxX - utmMinX) / TD.Header.SizeX;
            double cellsizeY = (utmMaxY - utmMinY) / TD.Header.SizeY;

            newTD.Header.CellSizeX = cellsizeX;
            newTD.Header.CellSizeY = cellsizeY;

            //======================以上只是为了求出新文件的网格间距============================

            int sizeX = (int)Math.Round((maxx - minx) / cellsizeX);
            int sizeY = (int)Math.Round((maxy - miny) / cellsizeY);

            newTD.Header.SizeX = sizeX;
            newTD.Header.SizeY = sizeY;

            //====================以上是新文件的网格数=============================================

            //新文件的映射
            double[] geo = new double[6];
            geo[0] = minx;
            geo[1] = cellsizeX;
            geo[2] = 0;
            geo[3] = maxy;
            geo[4] = 0;
            geo[5] = 0 - Math.Abs(cellsizeY);

            newTD.Header.Affine = geo;

            //求新文件的高程数据
            int[] Data = new int[sizeX * sizeY];

            for (int j = 0; j < sizeY; j++)
            {
                for (int i = 0; i < sizeX; i++)
                {
                    double X = minx + i * cellsizeX;
                    double Y = miny + j * cellsizeY;
                    //转换成经纬度
                    lccToLonglat.TransformPoint(coorOut, X, Y, 0);

                    //插值算法，求出该点的高程值
                    Data[j * sizeX + i] = (int)GetPtElevation(TD, coorOut[0], coorOut[1]);
                }
            }

            newTD.Data2D = new dynamic[newTD.Header.SizeY, newTD.Header.SizeX];

            for (int i = 0; i < newTD.Header.SizeY; i++)
            {
                for (int j = 0; j < newTD.Header.SizeX; j++)
                {
                    newTD.Data2D[i,j] = Data[i * newTD.Header.SizeX + j];
                }
            }


            return newTD;
        }

        /// <summary>
        /// 根据经纬度范围从经纬度文件裁剪,其实这个不分投影的
        /// </summary>
        /// <param name="TD"></param>
        /// <param name="utmMinX"></param>
        /// <param name="utmMinY"></param>
        /// <param name="utmMaxX"></param>
        /// <param name="utmMaxY"></param>
        /// <returns></returns>
        public static TerrainData Clip(TerrainData TD, double utmMinX, double utmMinY, double utmMaxX, double utmMaxY)
        {
            if (utmMinX >= TD.Header.CoordinateNE.X || utmMinY >= TD.Header.CoordinateNE.Y || utmMaxX <= TD.Header.CoordinateSW.X
                || utmMaxY <= TD.Header.CoordinateSW.Y)
            {
                //MessageBox.Show("网格不在地形文件的范围内！");
                return null;
            }

            if (utmMinX >= utmMaxX || utmMinY >= utmMaxY)
            {
                MessageBox.Show("坐标设置反了！");
                return null;
            }



            TerrainData newTD = new TerrainData();

            //新文件的范围坐标
            double minx = 0, miny = 0, maxx = 0, maxy = 0;

            minx = Math.Max(utmMinX, TD.Header.CoordinateSW.X);
            miny = Math.Max(utmMinY, TD.Header.CoordinateSW.Y);
            maxx = Math.Min(utmMaxX, TD.Header.CoordinateNE.X);
            maxy = Math.Min(utmMaxY, TD.Header.CoordinateNE.Y);

            newTD.Header.CoordinateSW = new HAPoint(minx, miny, 0);
            newTD.Header.CoordinateNW = new HAPoint(minx, maxy, 0);
            newTD.Header.CoordinateNE = new HAPoint(maxx, maxy, 0);
            newTD.Header.CoordinateSE = new HAPoint(maxx, miny, 0);

            //新文件的投影
            newTD.Header.Projection = TD.Header.Projection;

            //新文件的网格间距
            newTD.Header.CellSizeX = TD.Header.CellSizeX;
            newTD.Header.CellSizeY = TD.Header.CellSizeY;

            //新文件的映射
            double[] geo = new double[6];
            geo[0] = minx;
            geo[1] = newTD.Header.CellSizeX;
            geo[2] = 0;
            geo[3] = maxy;
            geo[4] = 0;
            geo[5] = 0 - Math.Abs(newTD.Header.CellSizeY);

            newTD.Header.Affine = geo;

            //无效数据
            newTD.Header.NoDataValue = TD.Header.NoDataValue;

            //新文件的网格数
            int startX = 0, startY = 0, endX = 0, endY = 0;
            startX = (int)Math.Floor((minx - TD.Header.CoordinateSW.X) / TD.Header.CellSizeX);
            startY = (int)Math.Floor((TD.Header.CoordinateNW.Y - maxy) / TD.Header.CellSizeY);
            endX = (int)Math.Floor((maxx - TD.Header.CoordinateSW.X) / TD.Header.CellSizeX) + 1;
            endY = (int)Math.Floor((TD.Header.CoordinateNW.Y - miny) / TD.Header.CellSizeY) + 1;

            newTD.Header.SizeX = endX - startX;
            newTD.Header.SizeY = endY - startY;

            int[] Data = new int[newTD.Header.SizeX * newTD.Header.SizeY];

            for (int j = 0; j < newTD.Header.SizeY; j++)
            {
                for (int i = 0; i < newTD.Header.SizeX; i++)
                {
                    Data[j * newTD.Header.SizeX + i] = TD.Data2D[startY + j,startX + i];
                }
            }

            newTD.Data2D = new dynamic[newTD.Header.SizeY, newTD.Header.SizeX];

            for (int i = 0; i < newTD.Header.SizeY; i++)
            {
                for (int j = 0; j < newTD.Header.SizeX; j++)
                {
                    newTD.Data2D[i,j] = Data[i * newTD.Header.SizeX + j];
                }
            }


            return newTD;

        }




        /// <summary>
        /// 经纬度文件转utm投影文件
        /// </summary>
        /// <param name="TD">地形文件数据</param>
        /// <param name="utmZone">utm区</param>
        /// <returns></returns>
        public static TerrainData ReprojectToUTM(TerrainData TD, int utmZone)
        {
            TerrainData newData = new TerrainData();

            //注意：在这里用到了两个坐标转换系统，一个是GDAL的，一个是Mapwindow的。
            //GDAL的速度快，所以首选是GDAL的。但是GDAL不支持经纬度坐标转特定区的UTM坐标，所以用到了Mapwindow的转换方法


            //=======================================================================================================
            //定义一个UTM空间投影
            OSGeo.OSR.SpatialReference utmPrj = new OSGeo.OSR.SpatialReference("");
            utmPrj.SetProjCS("UTM " + utmZone.ToString() + "(WGS84)");
            utmPrj.SetWellKnownGeogCS("WGS84");
            if (TD.Header.CoordinateSW.Y > 0)
            {
                utmPrj.SetUTM(utmZone, 1); //北半球
            }
            else
            {
                utmPrj.SetUTM(utmZone, 0); //南半球
            }
            //经纬度空间投影
            OSGeo.OSR.SpatialReference geoPrj = new OSGeo.OSR.SpatialReference(TD.Header.Projection);

            //坐标转换工具
            OSGeo.OSR.CoordinateTransformation utmToLonglat = new OSGeo.OSR.CoordinateTransformation(utmPrj, geoPrj);
            OSGeo.OSR.CoordinateTransformation longlatToUtm = new OSGeo.OSR.CoordinateTransformation(geoPrj, utmPrj);

            //以上是GDAL自待的转换方法
            //========================================================================================================
            //新文件的投影
            string prj = "";
            utmPrj.ExportToWkt(out prj);
            newData.Header.Projection = prj;

            //======================================================================================================


            double utmx, utmy;
            //int utmZone = 0;

            utmx = TD.Header.CoordinateSW.X;
            utmy = TD.Header.CoordinateSW.Y;

            double[] CoorOut = new double[3];

            //经纬度转换为utm坐标
            //西南角坐标
            //gdal方法
            longlatToUtm.TransformPoint(CoorOut, utmx, utmy, 0);
            HAPoint ptSW = new HAPoint(CoorOut[0], CoorOut[1], 0);
            //map6方法
            //HaGisPack.CoordinateConvert.LongLatToUTM(ref utmx, ref utmy, utmZone);
            //HAPoint ptSW = new HAPoint(utmx, utmy, 0);


            //西北角坐标
            utmx = TD.Header.CoordinateNW.X;
            utmy = TD.Header.CoordinateNW.Y;

            longlatToUtm.TransformPoint(CoorOut, utmx, utmy, 0);
            HAPoint ptNW = new HAPoint(CoorOut[0], CoorOut[1], 0);

            //HaGisPack.CoordinateConvert.LongLatToUTM(ref utmx, ref utmy, utmZone);
            //HAPoint ptNW = new HAPoint(utmx, utmy, 0);


            //东北角坐标
            utmx = TD.Header.CoordinateNE.X;
            utmy = TD.Header.CoordinateNE.Y;

            longlatToUtm.TransformPoint(CoorOut, utmx, utmy, 0);
            HAPoint ptNE = new HAPoint(CoorOut[0], CoorOut[1], 0);

            //HaGisPack.CoordinateConvert.LongLatToUTM(ref utmx, ref utmy, utmZone);
            //HAPoint ptNE = new HAPoint(utmx, utmy, 0);


            //东南角坐标
            utmx = TD.Header.CoordinateSE.X;
            utmy = TD.Header.CoordinateSE.Y;

            longlatToUtm.TransformPoint(CoorOut, utmx, utmy, 0);
            HAPoint ptSE = new HAPoint(CoorOut[0], CoorOut[1], 0);

            //HaGisPack.CoordinateConvert.LongLatToUTM(ref utmx, ref utmy, utmZone);
            //HAPoint ptSE = new HAPoint(utmx, utmy, 0);

            double utmMinX = 0, utmMaxX = 0, utmMinY = 0, utmMaxY = 0;

            utmMinX = Math.Min(ptNW.X, ptSW.X);
            utmMaxX = Math.Max(ptNE.X, ptSE.X);
            utmMinY = Math.Min(ptSW.Y, ptSE.Y);
            utmMaxY = Math.Max(ptNW.Y, ptNE.Y);

            //四个点坐标
            newData.Header.CoordinateSW = new HAPoint(utmMinX, utmMinY, 0);
            newData.Header.CoordinateNW = new HAPoint(utmMinX, utmMaxY, 0);
            newData.Header.CoordinateNE = new HAPoint(utmMaxX, utmMaxY, 0);
            newData.Header.CoordinateSE = new HAPoint(utmMaxX, utmMinY, 0);


            //===============原文件的坐标转换完毕==========================================


            //求新文件的网格间距与网格数
            double cellsizeX = (utmMaxX - utmMinX) / TD.Header.SizeX;
            double cellsizeY = (utmMaxY - utmMinY) / TD.Header.SizeY;

            newData.Header.CellSizeX = cellsizeX;
            newData.Header.CellSizeY = cellsizeY;

            newData.Header.SizeX = TD.Header.SizeX;
            newData.Header.SizeY = TD.Header.SizeY;

            //==============================================================================

            //新文件的映射
            double[] geo = new double[6];
            geo[0] = utmMinX;
            geo[1] = cellsizeX;
            geo[2] = 0;
            geo[3] = utmMaxY;
            geo[4] = 0;
            geo[5] = 0 - Math.Abs(cellsizeY);

            newData.Header.Affine = geo;

            //新文件的无效数据
            newData.Header.NoDataValue = TD.Header.NoDataValue;

            //求新文件的高程数据
            int[] Data = new int[TD.Header.SizeX * TD.Header.SizeY];

            double[] coor = new double[3];

            for (int j = 0; j < TD.Header.SizeY; j++)
            {
                for (int i = 0; i < TD.Header.SizeX; i++)
                {
                    double X = utmMinX + i * cellsizeX;
                    double Y = utmMaxY - j * cellsizeY;
                    //转换成经纬度
                    //HaGisPack.CoordinateConvert.UTMToLongLat(ref X, ref Y, utmZone, true);
                    //两种转换的方法，下面这种更快
                    utmToLonglat.TransformPoint(coor, X, Y, 0);
                    X = coor[0];
                    Y = coor[1];

                    //由于转换坐标图形会发生变化，再转换后的经纬度可能不在原来的范围内
                    if (X < TD.Header.CoordinateSW.X || X > TD.Header.CoordinateNE.X || Y < TD.Header.CoordinateSW.Y || Y > TD.Header.CoordinateNE.Y)
                    {
                        Data[j * TD.Header.SizeX + i] = (int)TD.Header.NoDataValue;
                    }
                    else
                    {
                        //插值算法，求出该点的高程值
                        Data[j * TD.Header.SizeX + i] = (int)GetPtElevation(TD, X, Y);
                    }
                }
            }

            newData.Data2D = new dynamic[newData.Header.SizeY, newData.Header.SizeX];

            for (int i = 0; i < newData.Header.SizeY; i++)
            {
                for (int j = 0; j < newData.Header.SizeX; j++)
                {
                    newData.Data2D[i,j] = Data[i * newData.Header.SizeX + j];
                }
            }


            return newData;
        }

        /// <summary>
        /// 经纬度文件转lcc投影文件
        /// </summary>
        /// <param name="TD">地形文件数据</param>
        /// <param name="lat_1">标准纬线1</param>
        /// <param name="lat_2">标准纬线2</param>
        /// <param name="originLong">中央经线</param>
        /// <param name="origionLat">（好像是）中央纬线</param>
        /// <param name="falseEasting">东偏移量</param>
        /// <param name="falseNorthing">西偏移量</param>
        /// <returns></returns>
        public static TerrainData ReprojectToLCC(TerrainData TD, double lat_1, double lat_2,
            double originLong, double origionLat, double falseEasting, double falseNorthing)
        {
            TerrainData newData = new TerrainData();

            //定义一个兰伯特投影
            OSGeo.OSR.SpatialReference lccPrj = new OSGeo.OSR.SpatialReference("");
            lccPrj.SetProjCS("LCC (WGS84)");
            lccPrj.SetWellKnownGeogCS("WGS84");
            lccPrj.SetLCC(lat_1, lat_2, origionLat, originLong, falseEasting, falseNorthing);


            //经纬度空间投影
            OSGeo.OSR.SpatialReference geoPrj = new OSGeo.OSR.SpatialReference(TD.Header.Projection);

            //坐标转换工具
            OSGeo.OSR.CoordinateTransformation lccToLonglat = new OSGeo.OSR.CoordinateTransformation(lccPrj, geoPrj);
            OSGeo.OSR.CoordinateTransformation longlatToLcc = new OSGeo.OSR.CoordinateTransformation(geoPrj, lccPrj);

            //以上是GDAL自待的转换方法
            //========================================================================================================
            //新文件的投影
            string prj = "";
            lccPrj.ExportToWkt(out prj);
            newData.Header.Projection = prj;

            //======================================================================================================


            double utmx, utmy;
            //int utmZone = 0;

            utmx = TD.Header.CoordinateSW.X;
            utmy = TD.Header.CoordinateSW.Y;

            double[] CoorOut = new double[3];

            //经纬度转换为utm坐标
            //西南角坐标
            //gdal方法
            longlatToLcc.TransformPoint(CoorOut, utmx, utmy, 0);
            HAPoint ptSW = new HAPoint(CoorOut[0], CoorOut[1], 0);
            //map6方法
            //HaGisPack.CoordinateConvert.LongLatToUTM(ref utmx, ref utmy, utmZone);
            //HAPoint ptSW = new HAPoint(utmx, utmy, 0);


            //西北角坐标
            utmx = TD.Header.CoordinateNW.X;
            utmy = TD.Header.CoordinateNW.Y;

            longlatToLcc.TransformPoint(CoorOut, utmx, utmy, 0);
            HAPoint ptNW = new HAPoint(CoorOut[0], CoorOut[1], 0);

            //HaGisPack.CoordinateConvert.LongLatToUTM(ref utmx, ref utmy, utmZone);
            //HAPoint ptNW = new HAPoint(utmx, utmy, 0);


            //东北角坐标
            utmx = TD.Header.CoordinateNE.X;
            utmy = TD.Header.CoordinateNE.Y;

            longlatToLcc.TransformPoint(CoorOut, utmx, utmy, 0);
            HAPoint ptNE = new HAPoint(CoorOut[0], CoorOut[1], 0);

            //HaGisPack.CoordinateConvert.LongLatToUTM(ref utmx, ref utmy, utmZone);
            //HAPoint ptNE = new HAPoint(utmx, utmy, 0);


            //东南角坐标
            utmx = TD.Header.CoordinateSE.X;
            utmy = TD.Header.CoordinateSE.Y;

            longlatToLcc.TransformPoint(CoorOut, utmx, utmy, 0);
            HAPoint ptSE = new HAPoint(CoorOut[0], CoorOut[1], 0);

            //HaGisPack.CoordinateConvert.LongLatToUTM(ref utmx, ref utmy, utmZone);
            //HAPoint ptSE = new HAPoint(utmx, utmy, 0);

            double utmMinX = 0, utmMaxX = 0, utmMinY = 0, utmMaxY = 0;

            utmMinX = Math.Min(ptNW.X, ptSW.X);
            utmMaxX = Math.Max(ptNE.X, ptSE.X);
            utmMinY = Math.Min(ptSW.Y, ptSE.Y);
            utmMaxY = Math.Max(ptNW.Y, ptNE.Y);

            //四个点坐标
            newData.Header.CoordinateSW = new HAPoint(utmMinX, utmMinY, 0);
            newData.Header.CoordinateNW = new HAPoint(utmMinX, utmMaxY, 0);
            newData.Header.CoordinateNE = new HAPoint(utmMaxX, utmMaxY, 0);
            newData.Header.CoordinateSE = new HAPoint(utmMaxX, utmMinY, 0);


            //===============原文件的坐标转换完毕==========================================


            //求新文件的网格间距与网格数
            double cellsizeX = (utmMaxX - utmMinX) / TD.Header.SizeX;
            double cellsizeY = (utmMaxY - utmMinY) / TD.Header.SizeY;

            newData.Header.CellSizeX = cellsizeX;
            newData.Header.CellSizeY = cellsizeY;

            newData.Header.SizeX = TD.Header.SizeX;
            newData.Header.SizeY = TD.Header.SizeY;

            //==============================================================================

            //新文件的映射
            double[] geo = new double[6];
            geo[0] = utmMinX;
            geo[1] = cellsizeX;
            geo[2] = 0;
            geo[3] = utmMaxY;
            geo[4] = 0;
            geo[5] = 0 - Math.Abs(cellsizeY);

            newData.Header.Affine = geo;

            //新文件的无效数据
            newData.Header.NoDataValue = TD.Header.NoDataValue;

            //求新文件的高程数据
            int[] Data = new int[TD.Header.SizeX * TD.Header.SizeY];

            double[] coor = new double[3];

            for (int j = 0; j < TD.Header.SizeY; j++)
            {
                for (int i = 0; i < TD.Header.SizeX; i++)
                {
                    double X = utmMinX + i * cellsizeX;
                    double Y = utmMinY + j * cellsizeY;
                    //转换成经纬度
                    //HaGisPack.CoordinateConvert.UTMToLongLat(ref X, ref Y, utmZone, true);
                    //两种转换的方法，下面这种更快
                    lccToLonglat.TransformPoint(coor, X, Y, 0);
                    X = coor[0];
                    Y = coor[1];

                    //由于转换坐标图形会发生变化，再转换后的经纬度可能不在原来的范围内
                    if (X < TD.Header.CoordinateSW.X || X > TD.Header.CoordinateNE.X || Y < TD.Header.CoordinateSW.Y || Y > TD.Header.CoordinateNE.Y)
                    {
                        Data[j * TD.Header.SizeX + i] = (int)TD.Header.NoDataValue;
                    }
                    else
                    {
                        //插值算法，求出该点的高程值
                        Data[j * TD.Header.SizeX + i] = (int)GetPtElevation(TD, X, Y);
                    }
                }
            }

            newData.Data2D = new dynamic[newData.Header.SizeY, newData.Header.SizeX];

            for (int i = 0; i < newData.Header.SizeY; i++)
            {
                for (int j = 0; j < newData.Header.SizeX; j++)
                {
                    newData.Data2D[i,j] = Data[i * newData.Header.SizeX + j];
                }
            }

            return newData;

        }


        # endregion



        # region  工具算法



        /// <summary>
        /// 求任意一点的高程值
        /// </summary>
        /// <param name="TD"></param>
        /// <param name="X">经度</param>
        /// <param name="Y">纬度</param>
        /// <returns></returns>
        public static double GetPtElevation(TerrainData TD, double X, double Y)
        {
            double h = 0;

            //先判断坐标是否超出了本文件的坐标范围
            if (X < TD.Header.CoordinateSW.X || Y < TD.Header.CoordinateSW.Y || X > TD.Header.CoordinateNE.X ||
                Y > TD.Header.CoordinateNE.Y)
            {
                h = TD.Header.NoDataValue;
            }
            else
            {
                double startX = TD.Header.CoordinateNW.X, startY = TD.Header.CoordinateNW.Y;
                HAPoint ptsw = new HAPoint();
                HAPoint ptnw = new HAPoint();
                HAPoint ptne = new HAPoint();
                HAPoint ptse = new HAPoint();

                //找出邻近的四个点坐标

                int i = 0;
                int j = 0;


                i = (int)Math.Floor((X - TD.Header.CoordinateNW.X) / TD.Header.CellSizeX);
                j = (int)Math.Floor((TD.Header.CoordinateNW.Y - Y) / TD.Header.CellSizeY);

                //如果出界了
                if (i >= TD.Header.SizeX - 1)
                {
                    if (j < TD.Header.SizeY - 1)
                    {
                        h = TD.Data2D[j,TD.Header.SizeX - 1];
                    }
                    else
                    {
                        h = TD.Data2D[TD.Header.SizeY - 1,TD.Header.SizeX - 1];
                    }
                    return h;
                }
                else
                {
                    if (j >= TD.Header.SizeY - 1)
                    {
                        h = TD.Data2D[TD.Header.SizeY - 1,i];
                        return h;
                    }
                }


                //第二种方法
                ptsw.X = startX + i * TD.Header.CellSizeX;
                ptnw.X = startX + i * TD.Header.CellSizeX;
                ptne.X = startX + (i + 1) * TD.Header.CellSizeX;
                ptse.X = startX + (i + 1) * TD.Header.CellSizeX;

                ptsw.Y = startY - (j + 1) * TD.Header.CellSizeY;
                ptse.Y = startY - (j + 1) * TD.Header.CellSizeY;
                ptne.Y = startY - j * TD.Header.CellSizeY;
                ptnw.Y = startY - j * TD.Header.CellSizeY;

                //第一种求位置法:遍历法
                //for (i = 0; i < this.SizeX; i++)
                //{
                //    if (startX > X)
                //    {
                //        ptsw.X = startX - this.CellSizeX;
                //        ptnw.X = startX - this.CellSizeX;
                //        ptne.X = startX;
                //        ptse.X = startX;
                //        break;
                //    }
                //    startX += this.CellSizeX;
                //}

                //for (j = 0; j < this.SizeY; j++)
                //{
                //    if (startY > Y)
                //    {
                //        ptsw.Y = startY - this.CellSizeY;
                //        ptse.Y = startY - this.CellSizeY;
                //        ptne.Y = startY;
                //        ptnw.Y = startY;
                //        break;
                //    }
                //    startY += this.CellSizeY;
                //}

                //四个点的Z值
                //if (j < this.SizeY - 1)
                //{
                if ((j + 1) * TD.Header.SizeX + i < TD.Header.SizeX * TD.Header.SizeY)
                {
                    ptsw.Z = TD.Data2D[j + 1,i];
                }
                else
                {
                    ptsw.Z = 0;
                }


                if (j * TD.Header.SizeX + i < TD.Header.SizeX * TD.Header.SizeY)
                {
                    ptnw.Z = TD.Data2D[j,i];
                }
                else
                {
                    ptnw.Z = 0;
                }


                if (j * TD.Header.SizeX + i + 1 < TD.Header.SizeX * TD.Header.SizeY)
                {

                    ptne.Z = TD.Data2D[j,i + 1];
                }
                else
                {
                    ptne.Z = 0;
                }


                if ((j + 1) * TD.Header.SizeX + i + 1 < TD.Header.SizeX * TD.Header.SizeY)
                {
                    ptse.Z = TD.Data2D[j + 1,i + 1];
                }
                else
                {
                    ptse.Z = 0;
                }


                HAPoint[] pt = new HAPoint[4];

                pt[0] = ptsw;
                pt[1] = ptnw;
                pt[2] = ptne;
                pt[3] = ptse;

                if (ptsw.Z == TD.Header.NoDataValue && ptnw.Z == TD.Header.NoDataValue &&
                    ptne.Z == TD.Header.NoDataValue && ptse.Z == TD.Header.NoDataValue)
                {
                    h = TD.Header.NoDataValue;
                    return h;
                }

                if (ptsw.Z == TD.Header.NoDataValue) ptsw.Z = 0;
                if (ptnw.Z == TD.Header.NoDataValue) ptnw.Z = 0;
                if (ptne.Z == TD.Header.NoDataValue) ptne.Z = 0;
                if (ptse.Z == TD.Header.NoDataValue) ptse.Z = 0;

                //插值算法
                h = Interpolation(X, Y, pt);
                //}
                //else
                //{
                //    h = this.NoDataValue;
                //}
            }
            return h;

        }




        ///// <summary>
        ///// 求任意一点的高程值
        ///// </summary>
        ///// <param name="TD"></param>
        ///// <param name="X">经度</param>
        ///// <param name="Y">纬度</param>
        ///// <returns></returns>
        //public static double GetPtElevation(TerrainData TD, double X, double Y)
        //{
        //    double h = 0;

        //    //先判断坐标是否超出了本文件的坐标范围
        //    if (X < TD.Header.CoordinateSW.X || Y < TD.Header.CoordinateSW.Y || X > TD.Header.CoordinateNE.X ||
        //        Y > TD.Header.CoordinateNE.Y)
        //    {
        //        h = TD.Header.NoDataValue;
        //    }
        //    else
        //    {
        //        double startX = TD.Header.CoordinateNW.X, startY = TD.Header.CoordinateNW.Y;
        //        HAPoint ptsw = new HAPoint();
        //        HAPoint ptnw = new HAPoint();
        //        HAPoint ptne = new HAPoint();
        //        HAPoint ptse = new HAPoint();

        //        //找出邻近的四个点坐标

        //        int i = 0;
        //        int j = 0;


        //        i = (int)Math.Floor((X - TD.Header.CoordinateNW.X) / TD.Header.CellSizeX);
        //        j = (int)Math.Floor((TD.Header.CoordinateNW.Y - Y) / TD.Header.CellSizeY);

        //        //第二种方法
        //        ptsw.X = startX + i * TD.Header.CellSizeX;
        //        ptnw.X = startX + i * TD.Header.CellSizeX;
        //        ptne.X = startX + (i + 1) * TD.Header.CellSizeX;
        //        ptse.X = startX + (i + 1) * TD.Header.CellSizeX;

        //        ptsw.Y = startY - (j + 1) * TD.Header.CellSizeY;
        //        ptse.Y = startY - (j + 1) * TD.Header.CellSizeY;
        //        ptne.Y = startY - j * TD.Header.CellSizeY;
        //        ptnw.Y = startY - j * TD.Header.CellSizeY;

        //        //第一种求位置法:遍历法
        //        //for (i = 0; i < this.SizeX; i++)
        //        //{
        //        //    if (startX > X)
        //        //    {
        //        //        ptsw.X = startX - this.CellSizeX;
        //        //        ptnw.X = startX - this.CellSizeX;
        //        //        ptne.X = startX;
        //        //        ptse.X = startX;
        //        //        break;
        //        //    }
        //        //    startX += this.CellSizeX;
        //        //}

        //        //for (j = 0; j < this.SizeY; j++)
        //        //{
        //        //    if (startY > Y)
        //        //    {
        //        //        ptsw.Y = startY - this.CellSizeY;
        //        //        ptse.Y = startY - this.CellSizeY;
        //        //        ptne.Y = startY;
        //        //        ptnw.Y = startY;
        //        //        break;
        //        //    }
        //        //    startY += this.CellSizeY;
        //        //}

        //        //四个点的Z值
        //        //if (j < this.SizeY - 1)
        //        //{
        //        if ((j + 1) * TD.Header.SizeX + i < TD.Header.SizeX * TD.Header.SizeY)
        //        {
        //            ptsw.Z = TD.Data2D[j + 1][i];
        //        }
        //        else
        //        {
        //            ptsw.Z = 0;
        //        }
        //        if (ptsw.Z == TD.Header.NoDataValue) ptsw.Z = 0;

        //        if (j * TD.Header.SizeX + i < TD.Header.SizeX * TD.Header.SizeY)
        //        {
        //            ptnw.Z = TD.Data2D[j][i];
        //        }
        //        else
        //        {
        //            ptnw.Z = 0;
        //        }
        //        if (ptnw.Z == TD.Header.NoDataValue) ptnw.Z = 0;

        //        if (j * TD.Header.SizeX + i + 1 < TD.Header.SizeX * TD.Header.SizeY)
        //        {
        //            ptne.Z = TD.Data2D[j][i + 1];
        //        }
        //        else
        //        {
        //            ptne.Z = 0;
        //        }
        //        if (ptne.Z == TD.Header.NoDataValue) ptne.Z = 0;

        //        if ((j + 1) * TD.Header.SizeX + i + 1 < TD.Header.SizeX * TD.Header.SizeY)
        //        {
        //            ptse.Z = TD.Data2D[j + 1][i + 1];
        //        }
        //        else
        //        {
        //            ptse.Z = 0;
        //        }
        //        if (ptse.Z == TD.Header.NoDataValue) ptse.Z = 0;

        //        HAPoint[] pt = new HAPoint[4];

        //        pt[0] = ptsw;
        //        pt[1] = ptnw;
        //        pt[2] = ptne;
        //        pt[3] = ptse;

        //        //插值算法
        //        h = Interpolation(X, Y, pt);
        //        //}
        //        //else
        //        //{
        //        //    h = this.NoDataValue;
        //        //}
        //    }
        //    return h;

        //}

        /// <summary>
        /// 插值算法
        /// </summary>
        /// <param name="x"></param>
        /// <param name="y"></param>
        /// <param name="pt"></param>
        /// <returns></returns>
        public static double Interpolation(Double x, Double y, HAPoint[] pt)
        {
            Double z1, z2, z3, z4;
            Double S = pt[1].Y - pt[0].Y;
            Double SZ = pt[1].Z - pt[0].Z;

            z1 = pt[0].Z + (y - pt[0].Y) * SZ / S;

            SZ = pt[2].Z - pt[3].Z;
            z2 = pt[3].Z + (y - pt[3].Y) * SZ / S;

            S = pt[2].X - pt[1].X;
            SZ = pt[2].Z - pt[1].Z;
            z3 = pt[1].Z + (x - pt[1].X) * SZ / S;


            SZ = pt[3].Z - pt[0].Z;
            z4 = pt[0].Z + (x - pt[0].X) * SZ / S;

            Double Z = (z1 + z2 + z3 + z4) / 4;
            return Z;


        }
        


        # endregion

    }
}
